Curso II – SISSA

Evaluación de Productos de Precipitación

Objetivo y contenido

Objetivo y contenido

Objetivo: explorar la importancia de evaluar productos de precipitación, consideraciones clave, funciones objetivo utilizadas y el preprocesamiento de datos necesario antes de la evaluación.

  • ¿Por qué evaluar productos de precipitación?
  • Consideraciones
  • Funciones objetivo
  • Preprocesamiento de datos
  • Ejemplo: evaluando productos de precipitación en el sur de África

¿Por qué evaluar productos de precipitación?

  • Los productos de precipitación son fundamentales para comprender y monitorear el ciclo del agua y su impacto en diversos sectores como la agricultura, hidrología, gestión de recursos hídricos, prevención de desastres, entre otros.

  • La evaluación de la calidad y precisión de estos productos es crucial para garantizar la confiabilidad de los datos y su utilidad en la toma de decisiones y aplicaciones relacionadas.

Consideraciones

Consideraciones

Resolución espacial y temporal (Benedict et al., 2019)

Consideraciones

El desempeño de los productos (Baez-Villanueva et al., 2020)

Consideraciones

¿Utilizarias el mismo producto para modelación hidrológica, cálculo de índices de sequía (SPI) y análisis de asignación de agua a largo plazo?

Consideraciones

Hora de reporte de las estaciones (Beck et al., 2019)

Evaluación de productos

Funciones objetivo

Las funciones objetivo son métricas utilizadas para evaluar el desempeño de los productos de precipitación. Algunas funciones objetivo comunes son:

  • RMSE
  • NRMSE
  • MAE
  • coeficiente de correlación (r)
  • NSE
  • KGE
  • PBIAS

Funciones objetivo

Adicionalmente, queremos analizar el desempeño de los productos en detectar ciertas intensidades de precipitación, para esto podemos usar funciones objetivo categoricas.

  • POD
  • FAR
  • fBIAS
  • CSI

Para una descripción mas detallada de estas funciones objetivo categóricas puedes leer Zambrano-Bigiarini et al., (2017) y Baez-Villanueva et al., (2018)

Funciones objetivo

En este módulo, utilizaremos el índice KGE’ (Kling-Gupta Efficiency) ya que es capaz de descomponer el rendimiento total en tres componentes distintos:

  • Correlación de Pearson
  • Bias ratio
  • Variability ratio

Funciones objetivo

\[ \textrm{KGE'} = 1 - \sqrt{ (r-1)^2 + (\beta-1)^2 + (\gamma-1)^2 } \]

\[ r = \frac{ \sum_{i=1}^n{(O_i - \bar{O}) (S_i - \bar{S}) } }{\sqrt{\sum_{i=1}^n{(O_i - \bar{O})^2}}\sqrt{\sum_{i=1}^n{(S_i - \bar{S})^2}}} \]

\[ \beta = \frac{\mu_{s}}{\mu_{o}} \]

\[ \gamma = \frac{CV_{s}}{CV_{o}} = \frac{\sigma_{s}/\mu_{s}}{\sigma_{o}/\mu_{o}} \]

Preprocesamiento de datos

Preprocesamiento de datos

Antes de la evaluación, es necesario realizar ciertos pasos de preprocesamiento en los datos de precipitación.

  • Por ejemplo: la limpieza y control de calidad de datos y la interpolación espacial y temporal para obtener una resolución uniforme.
  • El preprocesamiento adecuado garantiza que los datos estén en un formato consistente y listos para realizar la evaluación.

Ejemplo - Evaluacion del Desempeño de productos de precipitacion

Datos

Para este ejercicio, utilizaremos datos diarios de estaciones de la base de datos Global Historical Climatology Network daily (GHCNd) para el sur de África. La evaluación se realizará durante el periodo 1990-1992 y, para esto, también utilizaremos datos diarios de MSWEPv2.8 y CHIRPSv2.

Evaluación

Paso 1: Primero cargaremos las series temporales observadas y los metadatos. Para esto, almacenemos las rutas a los archivos:

observations_path <- 'C:/Users/Data/timeseries_1990-1992.zoo' 
metadata_path     <- 'C:/Users/Data/metadata.csv'

Después, utilizaremos la función read.zoo cargando el paquete hydroGOF para leer las series temporales, y la función read.csv para leer los metadatos.

library(hydroGOF)

Evaluación

observations <- read.zoo(observations_path, header = TRUE)
metadata     <- read.csv(metadata_path)
dim(observations)
## [1] 1096 1377
dim(metadata)
## [1] 1377    6

Podemos ver que tenemos un total de 1377 estaciones. Pero, ¿cuántos datos faltantes tienen estas estaciones?

Evaluación

observations[1:5, 1:5]
##            BC000068026 BC000068032 BC000068226 BC000068328 BC004612080
## 1990-01-01         0.0           0           0           0           0
## 1990-01-02         1.0           0           0           0           0
## 1990-01-03         0.5           0           0           0           0
## 1990-01-04        34.0           0           0           0           0
## 1990-01-05        14.5           0           0           0           0
head(metadata)
##            ID LATITUDE LONGITUDE ELEMENT YEAR_START YEAR_END
## 1 BC000068026  -18.367    21.850    PRCP       1932     1992
## 2 BC000068032  -19.983    23.417    PRCP       1921     2022
## 3 BC000068226  -24.017    21.883    PRCP       1958     1992
## 4 BC000068328  -26.050    22.450    PRCP       1939     1992
## 5 BC004612080  -26.467    20.617    PRCP       1960     2022
## 6 BC005452830  -25.220    25.670    PRCP       1921     1989

Evaluación

Paso 2: Analizar el número de datos faltantes en cada estación. Para esto, vamos a crear una función básica:

get_ratio_md <- function(x){
  
  md <- sum(is.na(x)) / length(x)
  return(md)
  
}

Evaluación

Ahora la aplicaremos a cada una de las series temporales (columnas del objeto zoo).

missing_data <- apply(observations, 2, get_ratio_md)
head(missing_data, 6)
##  BC000068026  BC000068032  BC000068226  BC000068328  BC004612080  BC005452830 
## 0.0027372263 0.0000000000 0.0009124088 0.0018248175 0.0000000000 1.0000000000

Evaluación

Podemos generar un histograma para visualizar la distribución de los datos faltantes:

hist(missing_data, breaks = 50)

Evaluación

Seleccionemos un umbral de 0.2 como límite de datos faltantes. Es decir, conservaremos las estaciones que tengan menos del 20% de datos faltantes.

exclude_pos      <- which(missing_data > 0.20)
exclude_stations <- names(observations)[exclude_pos]
length(exclude_stations)
## [1] 585
head(exclude_stations)
## [1] "BC005452830" "BC005836580" "BC006270240" "BC006282050" "BC007034890"
## [6] "BC007591160"

Evaluación

Ahora tenemos que remover dichas estaciones tanto de las series temporales como de los metadatos:

observations <- observations[,-exclude_pos]
dim(observations)
## [1] 1096  792

metadata_pos <- which(metadata$ID %in% exclude_stations)
metadata     <- metadata[-metadata_pos,]
dim(metadata)
## [1] 792   6

Evaluación

Paso 3: Una vez que tenemos las series temporales y los metadatos de las estaciones que utilizaremos, vamos a cargar los productos de precipitación. Para esto, guardemos las rutas de los productos en dos objetos:

chirps_path <- 'C:/Users/Data/CHIRPSv2/'
mswep_path  <- 'C:/Users/Data//MSWEPv2.8/'

Evaluación

Ahora listemos los archivos que se encuentran en ambas carpetas.

chirps <- list.files(chirps_path, full.names = TRUE)
mswep  <- list.files(mswep_path, full.names = TRUE)
head(chirps)
## [1] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900101.tif"
## [2] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900102.tif"
## [3] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900103.tif"
## [4] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900104.tif"
## [5] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900105.tif"
## [6] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900106.tif"
head(mswep)
## [1] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-01.tif"
## [2] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-02.tif"
## [3] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-03.tif"
## [4] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-04.tif"
## [5] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-05.tif"
## [6] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-06.tif"

Evaluación

Finalmente, podemos utilizar la función rast del paquete terra para importar cada conjunto de archivos:

chirps <- rast(chirps)
mswep <- rast(mswep)
print(chirps)
## class       : SpatRaster 
## dimensions  : 78, 86, 1096  (nrow, ncol, nlyr)
## resolution  : 0.25, 0.25  (x, y)
## extent      : 11.625, 33.125, -34.875, -15.375  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : chirps_daily_19900101.tif  
##               chirps_daily_19900102.tif  
##               chirps_daily_19900103.tif  
##               ... and 1093 more source(s)
## names       : chirp~00101, chirp~00102, chirp~00103, chirp~00104, chirp~00105, chirp~00106, ... 
## min values  :     0.00000,     0.00000,     0.00000,     0.00000,     0.00000,      0.0000, ... 
## max values  :    56.70765,    29.22468,    44.08159,    39.65398,    38.00595,     37.1624, ...
print(mswep)
## class       : SpatRaster 
## dimensions  : 78, 86, 1096  (nrow, ncol, nlyr)
## resolution  : 0.25, 0.25  (x, y)
## extent      : 11.625, 33.125, -34.875, -15.375  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : MSWEPv280_1990-01-01.tif  
##               MSWEPv280_1990-01-02.tif  
##               MSWEPv280_1990-01-03.tif  
##               ... and 1093 more source(s)
## names       : MSWEP~01-01, MSWEP~01-02, MSWEP~01-03, MSWEP~01-04, MSWEP~01-05, MSWEP~01-06, ... 
## min values  :     0.00000,     0.00000,     0.00000,     0.00000,     0.00000,     0.00000, ... 
## max values  :    53.84358,    23.13651,    28.26126,    45.11976,    37.28848,    48.36964, ...

Evaluación

Podemos graficar la primera capa de cada set de datos:

plot(chirps[[1]], main = 'CHIRPSv2')

Evaluación

Podemos graficar la primera capa de cada set de datos:

plot(mswep[[1]], main = 'MSWEPv2.8')

Evaluación

Paso 4: Convertir el archivo de metadatos en un objeto espacial para extraer la información de los productos en función de la ubicación de las estaciones.

head(metadata)
##            ID LATITUDE LONGITUDE ELEMENT YEAR_START YEAR_END
## 1 BC000068026  -18.367    21.850    PRCP       1932     1992
## 2 BC000068032  -19.983    23.417    PRCP       1921     2022
## 3 BC000068226  -24.017    21.883    PRCP       1958     1992
## 4 BC000068328  -26.050    22.450    PRCP       1939     1992
## 5 BC004612080  -26.467    20.617    PRCP       1960     2022
## 8 BC005847300  -24.667    25.917    PRCP       1925     1997
vct <- vect(metadata, geom=c("LONGITUDE", "LATITUDE"))
print(vct)
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 792, 4  (geometries, attributes)
##  extent      : 15.05, 32.42, -34.833, -17.817  (xmin, xmax, ymin, ymax)
##  coord. ref. :  
##  names       :          ID ELEMENT YEAR_START YEAR_END
##  type        :       <chr>   <chr>      <int>    <int>
##  values      : BC000068026    PRCP       1932     1992
##                BC000068032    PRCP       1921     2022
##                BC000068226    PRCP       1958     1992

Evaluación

plot(mswep[[1]], main = 'MSWEPv2.8')
points(vct)

Evaluación

Podemos usar el paquete maps para agregar el contorno de los paises:

library(maps)
plot(mswep[[1]], main = 'MSWEPv2.8')
points(vct)
map(add = TRUE)

Evaluación

Evaluación

Paso 5: Extraeremos la información correspondiente a la ubicación de las estaciones utilizando la función extract del paquete terra. Esta función nos permitirá extraer los valores de los productos de acuerdo a la ubicación de las estaciones espaciales previamente convertidas en un objeto espacial.

chirps_ts <- t(terra::extract(chirps, vct))
mswep_ts <- t(terra::extract(mswep, vct))
dim(chirps_ts)
## [1] 1097  792
dim(mswep_ts)
## [1] 1097  792

Evaluación

chirps_ts[1:5, 1:5]
##                            [,1]     [,2]         [,3] [,4] [,5]
## ID                    1.0000000 2.000000 3.0000000000    4    5
## chirps_daily_19900101 0.5458916 1.354978 0.0005568612    0    0
## chirps_daily_19900102 0.4041451 0.000000 0.0000000000    0    0
## chirps_daily_19900103 8.3321171 2.423540 0.0000000000    0    0
## chirps_daily_19900104 1.8144772 0.000000 0.0000000000    0    0
mswep_ts[1:5, 1:5]
##                            [,1]        [,2]        [,3]      [,4]       [,5]
## ID                    1.0000000 2.000000000 3.000000000 4.0000000 5.00000000
## MSWEPv280_1990-01-01  3.6173794 0.133561790 0.000000000 0.0000000 0.00000000
## MSWEPv280_1990-01-02  0.4801421 0.008661872 0.001714454 0.0000000 0.00000000
## MSWEPv280_1990-01-03  0.8579947 0.005163515 0.000000000 0.0000000 0.00000000
## MSWEPv280_1990-01-04 29.3101654 0.411440879 0.000000000 0.1994748 0.02766745

Evaluación

Removemos la primera fila de resultados en ambos data frames.

chirps_ts <- chirps_ts[-1,]
mswep_ts <- mswep_ts[-1,]
dim(chirps_ts)
## [1] 1096  792
dim(mswep_ts)
## [1] 1096  792

Evaluación

Podemos agregar el nombre de las estaciones a ambos data frames y convertirlos a objetos zoo.

colnames(chirps_ts) <- vct$ID
colnames(mswep_ts) <- vct$ID

chirps_ts <- zoo(chirps_ts, index(observations))
mswep_ts <- zoo(mswep_ts, index(observations))

Evaluación

Es importante destacar que la posición de las estaciones en los metadatos coincide con el orden de las extracciones espaciales. Es decir, la estación en la primera fila de los metadatos corresponde a la primera columna de las extracciones. Es fundamental asegurarse de que el orden de los metadatos sea consistente con el orden de las series temporales de las estaciones. En caso contrario, será necesario ordenarlos de manera correspondiente para garantizar la correcta correspondencia entre los metadatos y las extracciones.

Evaluación

Comparemos el orden de las estaciones entre metadatos y observaciones:

all(metadata$ID == names(observations))
## [1] TRUE

Evaluación

Para evaluar el desempeño de las estaciones, podemos utilizar la función KGE del paquete hydroGOF, como se muestra en el siguiente ejemplo:

KGE(sim = chirps_ts[,1], obs = observations[,1], method = '2012')
## [1] 0.3300503
KGE(sim = chirps_ts[,1], obs = observations[,1], method = '2012', out.type = 'full')
## $KGE.value
## [1] 0.3300503
## 
## $KGE.elements
##         r      Beta     Gamma 
## 0.4103995 0.9721374 0.6830971

Evaluación

Paso 5: En seguida, debemos de desarrollar un loop para evaluar cada uno de los productos con el KGE y sus componentes:

  • El ciclo tiene que iterar el numero total de estaciones analizadas
  • Debe de guardar el valor del KGE’ y sus componentes.

Evaluación

# Vectores vacios para almacenar los resultados
kges_chirps <- r_chirps <- beta_chirps <- gamma_chirps <- c()

for(i in 1:ncol(chirps_ts)){
  
  # Extraer los valores del producto y observados
  x <- observations[,i]
  y <- chirps_ts[,i]
 
  # Calcular el KGE' y sus componentes
  kge_all <- hydroGOF::KGE(sim = y, obs = x, method = "2012", out.type = "full")
  
  kges_chirps[i]     <- kge_all$KGE.value
  r_chirps[i]        <- kge_all$KGE.elements[1]
  beta_chirps[i]     <- kge_all$KGE.elements[2]
  gamma_chirps[i]    <- kge_all$KGE.elements[3]
  
}

# Generar un data frame con los resultados
  res_chirps <- data.frame(ID = names(observations), long = metadata$LONGITUDE, lat = metadata$LATITUDE, 
                    KGE = kges_chirps, r = r_chirps, beta = beta_chirps, gamma = gamma_chirps)

Evaluación

Podemos visualizar los primeros 10 elementos de los resultados:

head(res_chirps, 10)
##             ID   long     lat       KGE         r      beta     gamma
## 1  BC000068026 21.850 -18.367 0.3300503 0.4103995 0.9721374 0.6830971
## 2  BC000068032 23.417 -19.983 0.2477595 0.3962453 0.8472015 0.5781011
## 3  BC000068226 21.883 -24.017 0.3689620 0.4676964 0.9213105 0.6703491
## 4  BC000068328 22.450 -26.050 0.3453056 0.4882832 1.0760378 0.5987656
## 5  BC004612080 20.617 -26.467 0.2843022 0.4233796 1.1143524 0.5917650
## 6  BC005847300 25.917 -24.667 0.2703927 0.3324382 0.7880933 0.7955895
## 7  BC007137570 25.430 -23.120 0.3646196 0.4564381 0.7583529 0.7767167
## 8  BC008382520 21.650 -21.700 0.4369397 0.5456965 1.0085835 0.6674770
## 9  BC008948490 27.500 -21.217 0.4043519 0.5007794 0.8767525 0.6993582
## 10 BC012192590 25.150 -17.817 0.3782020 0.5599314 0.8304033 0.5947725

Ahora realizamos lo mismo para MSWEPv2.8.

Evaluación

# Vectores vacios para almacenar los resultados
kges_mswep <- r_mswep <- beta_mswep <- gamma_mswep <- c()

for(i in 1:ncol(mswep_ts)){
  
  # Extraer los valores del producto y observados
  x <- observations[,i]
  y <- mswep_ts[,i]
 
  # Calcular el KGE' y sus componentes
  kge_all <- hydroGOF::KGE(sim = y, obs = x, method = "2012", out.type = "full")
  
  kges_mswep[i]     <- kge_all$KGE.value
  r_mswep[i]        <- kge_all$KGE.elements[1]
  beta_mswep[i]     <- kge_all$KGE.elements[2]
  gamma_mswep[i]    <- kge_all$KGE.elements[3]
  
}

# Generar un data frame con los resultados
  res_mswep <- data.frame(ID = names(observations), long = metadata$LONGITUDE, lat = metadata$LATITUDE, 
                    KGE = kges_mswep, r = r_mswep, beta = beta_mswep, gamma = gamma_mswep)

Evaluación

Podemos visualizar los primeros 10 elementos de los resultados:

head(res_mswep, warning=FALSE, 10)
##             ID   long     lat       KGE         r      beta     gamma
## 1  BC000068026 21.850 -18.367 0.7732305 0.8600063 1.0388965 0.8258930
## 2  BC000068032 23.417 -19.983 0.3683645 0.4052135 1.0246169 0.7888451
## 3  BC000068226 21.883 -24.017 0.7068569 0.7975389 1.0938342 0.8099013
## 4  BC000068328 22.450 -26.050 0.5019049 0.6471822 1.2117271 0.7193046
## 5  BC004612080 20.617 -26.467 0.4202343 0.5808006 1.2333419 0.6744977
## 6  BC005847300 25.917 -24.667 0.5994549 0.7160007 1.0794620 0.7289529
## 7  BC007137570 25.430 -23.120 0.7497016 0.9218014 0.7961388 0.8776282
## 8  BC008382520 21.650 -21.700 0.5602473 0.6129723 0.9501728 0.7972458
## 9  BC008948490 27.500 -21.217 0.8897841 0.9970764 0.9166046 0.9279986
## 10 BC012192590 25.150 -17.817 0.8900011 0.9704940 0.9561920 0.9035116

Evaluación

Finalmente podemos comparar el KGE’ y sus componentes:

kges <- data.frame(CHIRPSv2 = res_chirps$KGE, MSWEPv28 = res_mswep$KGE)
boxplot(kges, ylim = c(0, 1), 
        main = "Desempeño de acuerdo al KGE'", col = c('cyan', 'purple'), ylab = "KGE'")

Evaluación

r <- data.frame(CHIRPSv2 = res_chirps$r, MSWEPv28 = res_mswep$r)
boxplot(kges, ylim = c(0, 1), 
        main = "Desempeño de acuerdo al r", col = c('cyan', 'purple'), ylab = "r")

Evaluación

beta <- data.frame(CHIRPSv2 = res_chirps$beta, MSWEPv28 = res_mswep$beta)
boxplot(kges, ylim = c(-1, 1), 
        main = "Desempeño de acuerdo al Beta", col = c('cyan', 'purple'), ylab = expression(beta))

Evaluación

gamma <- data.frame(CHIRPSv2 = res_chirps$gamma, MSWEPv28 = res_mswep$gamma)
boxplot(kges, ylim = c(-1, 1), 
        main = "Desempeño de acuerdo al Gamma", col = c('cyan', 'purple'), ylab = expression(gamma))

¡Gracias por su atención!